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^ . Abstract. We revisit the accuracy of the variational Monte Carlo (VMC) method by taking 

C/3 ' an example of ground state properties for the one-dimensional Hubbard model. We start from 

the variational wave functions with the Gutzwiller and long-range Jastrow factor introduced 
J2 ■ by Capello et al. [Phys. Rev. B 72, 085121 (2005)] and further improve it by considering 

several quantum-number projections and a generalized one-body wave function. We find that 
the quantum spin projection and total momentum projection greatly improve the accuracy of 
the ground state energy within 0.5% error, for both small and large systems at half filling. 
Besides, the momentum distribution function n{k) at quarter filling calculated up to 196 sites 
allows us direct estimate of the critical exponents of the charge correlations from the power-law 
behavior of n{k) near the Fermi wave vector. Estimated critical exponents well reproduce those 
predicted by the Tomonaga-Luttinger theory. 



o , 

\^ , 1. Introduction 

Strongly correlated electron systems have attracted much interest over the past decades. 
. Although rigorous ground states are obtained analytically for some simple models, numerical 

I methods offer powerful tools for understanding more realistic systems. There are several 

numerical methods for the purpose of calculating low-energy properties of the strongly correlated 
systems. At finite temperatures, the auxiliary-field quantum Monte Carlo (AFQMC) method [H 
^ ' El E] gives exact ground states within the statistical errors. However, for the geometrically 
^ ■ frustrated systems, the errors become exponentially large, known as the negative sign problem. 

To overcome this difficulty, the Gaussian-basis Monte Carlo (GBMC) method [H El [6] has 
been proposed as an alternative to the AFQMC method, although the tractable range of 
parameters are still under debate. The dynamical mean-field theory (DMFT) [71 [8] can also 
be applied to strongly correlated electron systems. However, when we wish to take into account 
sufficient spatial correlations to represent fluctuations around broken symmetry states, we need 
larger clusters beyond the present computational accessibility. At zero temperature, the exact 
diagonalization (ED) method gives exact ground states. This method, however, is not able to 
handle large system sizes. The density matrix renormalization group (DMRG) method gives 
nearly exact ground states for relatively larger systems, though it is mainly applied to one- 
dimensional systems, or two-dimensional systems with cylindrical boundary conditions. The 
unbiased path-integral renormalization group (PIRG) method [TOl [HI [E] has been applied to 



several strongly correlated electron systems, with essentially accurate ground states, even in 
the presence of geometrical frustrations. However, strong correlation regime requires a set of 
many basis functions for accurate estimates, sometimes beyond the computational accessibility. 
Among them, although a bias inevitably exists, the VMC method [13] is applicable to large 
systems in the presence of strong correlations and geometrical frustrations. 

Accuracy of the calculated properties obtained by the VMC method relies on the choice 
of the trial wave functions. In most cases, the trial wave functions are given by Slater 
determinants \D) supplemented by correlation factors V as = V\D). One of the famous 
example of this trial function is a "projected Bardeen-Cooper-Schrieffer (BCS)" type wave 
function = "Pq" |BCS), where Vq = Y\i{^ ~ '^it'^^) denotes a Gutzwiller factor, which 
prohibits double occupancy of the sites, and |BCS) = nfc(''^fc + '^k<^k'\^~ki) I*-*) denotes a BCS 
state, which is the ground state of the BCS mean-field Hamiltonian. This projected BCS 
wave function is equivalent to the so-called "resonating valence bond (RVB)" wave function, 
which is the superposition of all dimer coverings \\.4\ IT^ . The projected BCS and RVB wave 
functions provide highly accurate descriptions in two-dimensional strongly correlated electron 
systems [121 E] . 

In principle, the trial wave functions based on Slater determinants become accurate for higher 
dimensions since the mean-field-like treatment becomes valid above the upper critical dimension. 
Although it seems difficult to capture the low-energy properties at lower spatial dimensions by 
using such a trial wave function, long-range correlation factors may restore the accuracy of the 
wave functions [18^ [T9] . Capello et al. have calculated the ground state of the one-dimensional 
Hubbard model by using the projected BCS type wave functions. They have shown that the 
long-range Jastrow factor is essential to describe the low-energy properties. They have also 
reproduced correct power-law behaviors of the spin and charge correlation functions. 

Although Capello et al.'s scheme appears to substantially reproduce the low-energy properties 
of the one-dimensional model, the errors in the energy become larger even for small system sizes, 
typically more than 2% for L = 18 for intermediate interaction strength [19]. This is mainly 
because they have used the Slater determinants with a limited number of variational parameters. 
The Slater determinants in their wave function are just a Fermi sea of the tight binding model, or 
short-range BCS wave functions with nonzero amplitudes at most up to third-nearest-neighbors. 
In addition, in practical numerical calculations, since the accessible system sizes are limited, it 
is important to prepare the wave functions preserve the correct quantum numbers within its 
accessible range of sizes. 

In this paper, to obtain more accurate wave functions of the low-dimensional systems for 
the smaller as well as for larger system sizes, we revisit the variational calculations of the one- 
dimensional Hubbard model, as a reference system. By considering the multi-variable Slater 
determinants and projections which restore the SU(2) spin rotational and spatially translational 
symmetries, we obtain more accurate ground states of the model than the previous studies. We 
also calculate the momentum distribution function n{k) up to 196 sites, the first high-precision 
variational calculation of n{k) to our knowledge, and directly estimate the critical exponents 
of the charge correlations. We compare them with the analysis by Kawakami et al. [2{)\ [21] as 
well as Ogata et al. |22j . Estimated critical exponents well reproduce those obtained by the 
Tomonaga-Luttinger theory. 

2. Model and method 

We consider the one-dimensional Hubbard model. The Hamiltonian is given by 




(1) 



where (cjo-), njo-, t and U denote the creation (annihilation) operator, the number operators, 
the nearest-neighbor hopping interaction, and the on-site Coulomb interaction, respectively. 

Exact ground states of the one-dimensional Hubbard model are obtained by Lieb and Wu [23|. 
At half filling, the ground state is a metal for U = 0, and is an insulator for nonzero positive U. 
Away from half filling, the ground state is a metal for any U >0. Unlike the normal Fermi liquid, 
the momentum distribution of this metal does not show a finite jump at the Fermi wave vector. 
Instead, it shows power-law decays near the Fermi wave vector. This type of the ground state, 
which appears in one-dimensional electron systems, is well known as the Tomonaga-Luttinger 
liquid [Ml EH]. 

To obtain accurate variational description of the ground state of the model, we use a multi- 
variable variational Monte Carlo (mVMC) method [T7]. We use the trial wave functions given 
by 

IV;) = VGVjVAhC^=^C''=^ IcApair) , (2) 

where \(t)pair), V,], Vdh, C^^'^, and C^^^ denote the one-body part, the Gutzwiller factor, 
the Jastrow factor, the doublon-holon factor, the spin quantum-number projection on to S" = 
space, and the total momentum projection on to -ftT = space, respectively, as we detail later. 
We calculate the ground states of the model under the periodic boundary conditions at half 
filling (n = 1) for N^^ = An + 2 sites, and at quarter filling [n = 0.5) for A'g = 8n + 4 sites, so 
that the systems obey the closed shell conditions. 

For the one-body part |i;^pair), we use a generalized pairing wave function defined as 



''pair/ 




(3) 



where fij, N^, and denote the variational parameters, the number of sites, and the number 
of electrons, respectively. When all /jj's are optimized simultaneously and independently, since 
the total number of fij is N^, the computation time for the optimization is highly demanding. 
To reduce the computation time to 0{Ns), we assume fij to have a sublattice structure, namely 
fij = fa{i){Ri ~ Rj), where a{i) denotes a sublattice index at site i. Since optimized ground 
states by using fij with reflectional symmetries, namely fij = fji, have higher energies, we do 
not impose reflectional symmetries on fij . Before calculating the ground states of larger system 
sizes, we examine the accuracy of the wave functions using several sublattice structures of fij for 
smaller system sizes less than about 20 sites. We find that optimized energies for two-sublattice 
structure and iVg-sublattice (fully optimized) structure are nearly the same. This suggests that 
the ground states are well-described by the wave functions with the two-sublattice structure for 
the one-dimensional Hubbard model. In the following, we show results for mainly two-sublattice 
structure, and one-sublattice translationally invariant structure. 

To include correlation effects beyond the one-body approximation, we consider the correlation 
factors 




Vg = exp \^-g2_^ni^niij , (4) 
Vj = exp ^Y^Vij{ni^ + nii){nj^ + nj^)^ , (5) 

and Vdh = exp - ^ ^ ami X] ' 



where g, Vij, and ami denote variational parameters of the Gutzwiher, Jastrow, and the doublon- 
holon correlation factors, respectively. We take into account long-range part of the Jastrow 
factor, which is indispensable for reproducing the low-energy properties of the one-dimensional 
Hubbard model |18tll9j. and the doublon-holon factor. A variable S,imi in Eq. ([6]) is a many-body 
operator, which gives one if a doublon (holon) exists at the i-th site and m holons (doublons) 
exist at l-th nearest neighbors, and gives zero otherwise. We assume the Jastrow factor to have 
translational symmetry, namely Vij = v{rij), where r^j = \Ri — Rj\ is the distance between 
the sites i and j. We note that the Gutzwiller factor and the doublon-holon factor have the 
translational symmetry by definition. 

In addition to the one-body part and the correlation factors, we introduce a projection C^^^, 
which restores the SU(2) spin rotational symmetry, and a projection C^^^, which restores 
the translational symmetry. The spin quantum- number and total momentum projections 
are switched on at all optimization steps. Since the correlation factors have the SU(2) 
and translational symmetries, the quantum-number projections commute with the correlation 
factors. 

To optimize a large number of these variational parameters, we use the stochastic 
reconfiguration (SR) method developed by Sorella [26]. We first optimize variational parameters 
fij for typically 1000 SR steps with fixed correlation factors. After confirming the energy 
convergence, we simultaneously optimize all the variational parameters (g, Vij, ami, and fij) 
for more than 4000 SR steps. 



3. Results 
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Figure 1. Calculated energy of Hubbard model at half filling. Statistical errors are smaller 
than the symbol sizes, (a) System size A's dependence at U/t = 4. The black crosses denote the 
exact energies obtained by the ED method for finite sizes and analytically at the thermodynamic 
limit. The purple down-pointing triangles and red squares denote calculated energies by the 
mVMC method with the translationally invariant and two-sublattice fij trial wave functions, 
respectively. The former is obtained without the quantum-number projections and the latter is 
with the spin and momentum projections to the total singlet and zero total momentum. Two 
and 0.5 percent errors measured from the exact results are plotted by dash-dotted purple and 
dashed red curves, respectively, (b) Coulomb interaction U/t dependence of errors in calculated 
energies at half filling for Ns = 18. The red squares denote our mVMC results obtained by using 
the two-sublattice fij trial wave functions. The errors in the energy are much smaller than their 
results. The green triangles and blue circles denote the errors in the energy obtained by Capello 
et al. by using short-range BCS wave functions with the third-nearest-neighbor amplitudes |19| . 
The former is obtained with the assumed Jastrow factor and the latter is with the optimized 
Jastrow factor. 



We obtain the ground-state energy of the one-dimensional Hubbard model at half filling, 
as shown in Fig. [TJ For the case of the two-sublattice fij, the errors in the energy are within 
approximately 0.5%, much smaller than the previous calculations by Capello et al. [19] The 
extended one-body part and quantum projections improve the accuracy of the ground state 
energy, even for the small system sizes. We show results of the two-sublattice fij condition, 
hereafter. 
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Figure 2. Calculated spin correlation functions of Hubbard model at U/t = 4, (a) at half 
filling and (b) at quarter filling, respectively. The red squares denote the mVMC results. These 
points are well described by the exact asymptotic expressions of the spin correlation functions, 
given by the dotted black lines. 



In the gapless one-dimensional systems, spin and charge correlations show power-law 
decays [27^ [28] . The asymptotic expressions of the spin correlation functions are given as 

C{r) = (5(0) • S{r)) = + Acos{2kFr)^^ + ■■■ , (7) 

where A, kp, and Kp denote a constant, the Fermi wave vector, and a critical exponent, 
respectively. The Fermi wave vector kp is written in terms of the filling n through the equation 
kp = n7r/2. The critical exponent Kp is a function of the Coulomb interaction U and the 
filling n. At half filling, the charge degrees of freedom are gapped out with the exponential 
decay of correlations in distance, leaving only the gapless spin correlation in the form of 
Eq. ([7]) formally with Kp = for nonzero repulsive U. Namely, the spin correlation functions 
decay as C(r) ~ A{—lY^/\nr /r \27\ \28\ 129] . Away from half filling, the critical exponent 
satisfies < Kp < 1 for nonzero repulsive U. Thus the spin correlation functions decay 
as Eq.Q. To include finite size effects under the periodic boundary conditions, we compare 
the calculated correlation functions with C(r) + C{Ns — r) instead of the original correlation 
function C(r) because of the refiection effects from the boundary. As shown in Fig. [21 calculated 
spin correlation functions of the Hubbard model reproduce \/ln r/r decay at half filling, and 
y/lnr /r^~^^p and decays at quarter filling with the critical exponent Kp ~ 0.7 at U/t = 4, 
respectively. This critical exponent is nearly the same as the exact value, Kp = 0.711, predicted 
by the Tomonaga-Luttinger theory [T8l [281 [20l [2T] . 

In our calculations, the one-body part of the trial wave function includes long-range part 
in fij and Vij. Although there should be no long-range order for the one-dimensional Hubbard 
model, our trial wave function could fallaciously stabilize the insulator with the antiferromagnetic 
long-range order. To examine whether our calculations correctly reproduce the absence of the 
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Figure 3. Peak values of spin structure factors of Hubbard model at \J jt = 4 and U/t = 10 
for each system size, (a) at half filling and (b) at quarter filling. Statistical errors are smaller 
than the symbol sizes. At half filling, the spin structure factor shows a peak at Qpeak = ^r, while 
at quarter filling, it shows a peak at Qpcak = 7r/2. The filled symbols denote results obtained by 
the ED method for the model under the periodic boundary conditions, while the open symbols 
denote the mVMC results. The red squares denote results at U/t = 4, and the green circles 
denote results at U/t = 10. The squared magnetizations, S{qpcak)/^s converge to zero in the 
thermodynamic limit. 

long-range order in the model, we calculate the peak values of the spin structure factors 

for each system size. For the insulator, size dependencies of the peak values of the spin structure 
factor are expected to follow the scaling 

Siqpc.^.) = ^J dr e^^p-^^ (S(0) • S{r)) ~ dr^ ~ (In N,f/\ (9) 
where A denotes the cutoff, while for the metal, since 

l/r^+Kp ^ i^^2 fQj. r^i^ 

we obtain 

5('?pcak) = lldr e'^p-^'^ (5(0) • 5(r)) ~ j^^ dr-j^ ~ const. + 0{N~'''). (10) 

By using these equations, we extrapolate the squared magnetization, the peak values of the 
spin structure factor «S'(q'peak) divided by the system size Ng, in the thermodynamic limit. As 
shown in Fig.[3l the squared magnetizations, S{qpcak)/]^s converge to zero in the thermodynamic 
limit. These results suggest that our trial wave functions with the long-range fij and symmetry- 
restoring projections C are able to reproduce ground states with no long-range magnetic order 
with power-law decay of correlations. 

For the one-dimensional Hubbard model, the momentum distribution 

-(^^) = ^E-^'^'''-''^-^(4c,.) (11) 

shows a power-law singularity 

\n{k) - n{kF)\ \k - kpl" (12) 




i u/t=m, n=o.5 

4' log A/7 = a log Mi+ iconst; 



\k - kpUn 

Figure 4. (a) Momentum distribution of Hubbard model at U/t = 10 for each system size, at 
quarter filling. Power-law singularities appear at the momenta kp = 7r/4 and Skp = 37r/4. (b) 
Logarithmic plot of momentum distribution near kp for each system size. Statistical errors are 
smaller than the symbol sizes. 



with a critical exponent a near the Fermi wave vector kp, contrary to the ordinary Fermi liquid, 
which shows a jump at the Fermi wave vector |201 121j . The critical exponent a is written in 
terms of Kp through the equation a = {Kp + 1/Kp — 2)/4 [30l [311 [28]. For the infinite Coulomb 
interactions, Ogata et al. have calculated the momentum distribution n{k) of the metal ground 
states by using the Bethe-ansatz wave functions, and have estimated the critical exponents a 
from Eq. (|12|) |22j . We calculate the momentum distribution of the Hubbard model at quarter 
filling up to 196 sites. As shown in Fig. Hl^a), calculated momentum distributions show power- 
law behavior near the wave vector k-p. Another singularity also appears at A; = 3A;F) as expected 
in the previous studies |22j . 
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Figure 5. Estimated critical exponent a from relation \n{k) — n{kp)\ oc \k — kp\°' near kp. 
Statistical errors are smaller than the symbol sizes. The red squares denote results at C//t = 4, 
and the green circles denote results at U/t = 10. These values are nearly the same as the exact 
results, a = 0.029 {U/t = 4) and a = 0.069 {U/t = 10), predicted by the Tomonaga-Luttinger 
theory [Ml [28l EOl EI] • 

To estimate the critical exponent a, we plot ln|?i(A;) — n(A;F)| as a function of \k — kp\ for 
each system size, as shown in Fig. Hl^b), and extrapolate the slopes from data points by using 



Eq. ()12p . Since large finite-size effects appear near \k — k-p\/TT < 0.1, we use data tliat satisfy 
\k — k-p\/TT G [0.1, 0.4] for fitting. As shown in Fig. [5l as tlie system size A'^g goes to infinity, the 
estimated critical exponent a(As) gradually converges to the exact value, a = 0.029 at U/t = A 
and a = 0.069 at U/t = 10, predicted by the Tomonaga-Luttinger theory [T8l [281 [201 [2T] . 



4. Summary 

To improve the VMC method, we revisit the variational description of the ground states of 
the one-dimensional Hubbard model, by considering the multi-variable Slater determinants as 
well as the long-range correlation factors. We find that the quantum spin projection and total 
momentum projection greatly improve the accuracy of the ground state energy with less than 
0.5% error, irrespective of the system size. By calculating the momentum distribution function 
n{k) of the model at quarter filling up to 196 sites, we have shown that our method is able to 
reproduce power-law decays at the momenta kp and 3kp. By directly estimating the critical 
exponents of the charge correlations from the power-law behavior of n(k) near kp, we have 
also shown that estimated critical exponents well reproduce those predicted by the Tomonaga- 
Luttinger theory. 
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